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ABSTRACT 


In this work, variations on velocity profiles in a flowing mass 
of molten glass in a forehearth are investigated. Formerly used 
parabolic velocity profiles are replaced with analytical solution of 
open channel flow equation, based on the available data on mass flow 
of molten glass through the channel in unit time. 

Concerning the viscosity effects; temperature dependence of 
viscosity is built in the model. However, it is assumed that the depth 
of the channel does not affect the viscosity gradients. 

To solve the system of non-linear Ferenc equations i.e., heat 
equation and flow equation; the analytical solution of the latter at 
the nodes is used for the numeric solution of the former iteratively, 
until the convergence is obtained. Predicted temperatures are compared 
to the available data from an actual operating forehearth, and against 
the results predicted by the previous model using simplying assumptions 


to prove its validity. 
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I. INTRODUCTION 


To improve the production of glass containers, and to transform the 
art of glass making into a science, the glass industry embarked an 
extensive forehearth measurement program in an attempt to resolve the 
complex phenomena of glass conditioning and cooling, while molten glass 
is taken out of the furnace and transferred to the forming machine. 

The objective of the program was the collection of quantitative data 
to achieve greater understanding of forehearth and to verify the math- 
ematical model of molten glass flow in forehearth developed by Duffin 
[Ref. 1]. 

A knowledge of the temperature distribution is needed for better under- 
Standing of the control process in glass plants. The above model proved 
that sophisticated models can give satisfactory results describing the 
forehearth behavior in a temperature sense. However they need a big 
computer and require considerable time to run. 

In this work, validity of the simplifying assumptions in [Ref. 1] 
are investigated and ways to relieve them are sought. Going to a more 
complicated model has academic interest to show that what was developed 


was a useful and fairly accurate model. 
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II. NATURE OF THE PROBLEM 


Forehearth is that section of the process where molten glass is 
transferred from furnace to the forming machine. During the flow, glass 
is conditioned in a temperature sense. A predetermined temperature is 
desired at the forming machine end, which assures minimum off-grade 
material after forming process. If undesirable temperature gradients 
exist at the input of forming machine, recycling of off-grade glassware 
is necessary, decreasing the overall system efficiency. 

What is needed is controlled conditioning throughout the forehearth 
while molten glass is flowing. This is mainly a heat transfer problem 
on a flowing mass of molten glass. 

To describe the physical phenomenon taking place in a typical fore- 
hearth used in glass container manufacture, see Figure (1). Glass flows 


from left to right in a open channel with dimensions approximately 26 


inches wide, 6 inches deep and 18 feet long. 
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Figure 1. Typical type K forehearth 
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Section AA is the point of entry from furnace area with average 
temperature being 2400 degrees F. at this point. Area from AA to BB 
is called cooling zone. Area from BB to CC is known as conditioning 
zone. This is the place where an attempt is made to condition the glass 
so that the delivery temperature would be constant in plane CC. After 
CC a glass gob is formed. A gob is a discrete mass of molten glass 
created by intermittenly shearing a stream of molten glass coming 
from an orifice. The gob formed falls into a guide chute, and is con- 
ducted to a forming machine mold. 

Temperature control is achieved by adjusting burner flame-levels and 
the amount of cooling wind in the cooling zone (section AA to BB); also 
with burner flame levels in the conditioning zone. 

Burners are mounted on manifolds and are located every 4.5 inches 
on both sides of the forehearth. They produce a blanket of heat over 
glass surface. Hot gases sweep over the ceramic radiating surface and 
heat exchange takes place between gases and surface which, in turn 
exchanges heat with molten glass. Cooling is achieved by introducing 
cooling air through inlet pipes. The wind exchanges heat with glass and 
the ceramic surface and goes out of the system through roof ports. Details 
of the forehearth cross section is shown in Figure (2). 

Heat exchange at the ceramic surface takes place by convection and 
radiation mechanisms. Heat is also conducted through the ceramic material 
in the sides and bottom of the forehearth. Possible disturbances can 
affect the system through these walls. 

Above discussion indicates a complex problem of heat exchange in fore- 
hearths, with boundaries as described. 

To describe temperature behavior within the glass, it is necessary to 


consider the mechanisms by which heat is transferred. 
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It had been shown that emission and absorption of radiation is a 
bulk phenomena. The interaction of simultaneous emission and absorption 
throughout the volume leads to radiative heat transfer between adjacent 
layers of material. When this radiative heat transfer is combined with 
true conductivity, mathematical description of this internal radiative 
mechanism becomes rather involved. 

Glass is, then, a diathermanous material whose diathermancy is 
determined by 


a) Wave length-absorption coefficient relation and 
b) Spectral composition of radiation involved. 


This dependency is usually shown in plots of absorption coefficient y 
versus wave length A, with temperature as a parameter. Any detailed 
analysis of the internal radiation should recognize this varying 
diathermancy. Integration can be replaced by summation using average 
values of y for various ranges of A, with the average value in a given 

A range being a function of temperature. Because of its transparency, 
glass layers beneath the surface will exchange heat with the surroundings. 
The depth to which this exchange occurs is appreciable varies with 

the type of the glass being considered. Upper one third of the glass flow- 
ing in a forehearth is considered as exchanging energy directly with the 
Surroundings by Duffin [Ref. 1] in his model. 

In the glass body, the radiant effects result from bulk phenomena of 
emission and absorption. Internally emitted radiation is reabsorbed by the 
glass directly and also as a result of internal reflections. This process 
leads to the "radiative conductance" of heat through the glass. Internal 
emission is characterized by “volume emissive power" which is the rate at 
which a unit volume of glass emits radiation in all directions. For an 


ideal gray material volume emissive power is given as, 
W = Ayn@oT" 


is 


Thus, it is proportional to absorption coefficient, refractive 
index and Stefan-Boltzman constant. According to Kellet [Ref. 2], in 
the interior of massive bodies of molten glass, steady-state temperature 
distributions turn out to be linear. It was also proposed that the 
effect of radiative conductivity be designated as an "radiation conduc- 
tivity". For the gray material it is given by | 


Bien st 


KRAD 3y 
Thus, the effective conductivity within the material is the sum of the 


thermal conductivity and radiation conductivity, as follows: 


k = k = k + 


eff true Kad 
This is set and held throughout the glass depth. But due to 
dependence on the absorption coefficient, variations are possible near 


the glass surface. 


A. DERIVATION OF DIFFERENTIAL EQUATIONS 

Derivation of the differential equations was accomplished by Duffin 
LRef. 1] as follows; using a right handed coordinate system and taking 
X coordinate as the flow direction, a differential volume element of 


dimensions dx, dy, dz is set up. Then writing energy balance with 


1. Energy Into Element Due To Mass Flow 
Area in X direction = dy-dz = dy 
Mass flow = V(1*dy)o 
Energy flow = (pdyV)(C)(T) 

2. Energy Out Of Element Due To Mass Flow 


e) 
ax (eCV, Tdy )dx ss oe CV Tdy 
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3. Net Flow Energy 
= Input - Output = 


,) 
- sx (eCV, Tdx)dy 


4, Energy Into Element Due To Conduction-Radiation 
Area = dx‘dz = dx 


oT 
ay 


5. Energy Out Of Element Due To Conduction-Radiation 


Energy = - k'dx — 


(- k'dx at) dy + [-k'dx at 


6. Net Conduction-Radiation Energy 
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7. Net Total Energy Flow 


9 = ee < a ae 
ay tk avi - sy (eo CV T)dxdy = (dxdy’p)C x 


so, differential equation is 
Joe: 


ay ae ax (CVT) = ec a5 


This being one dimensional model, actual case must include the wall 
effects which is not considered in above derivation. To account for that, 


a term for Z direction must be included. 


Kr (kt 2) - 2 (veer) = ec E 


Ihe - equal to zero, steady-state heat flow equation is obtained. 
In above equation density, specific heat and effective conductivity terms 
can be treated as temperature independent. Dependence of velocity on 


temperature will be discussed later. 
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B. BOUNDARY CONDITIONS 
To apply boundary conditions, a fixed coordinate system is placed 


in channel center as follows 





T @ x=0O y and z = 0 Temperature in the initial plane 

T @ yed x and z = 0 Radiating boundary 

T @ y=Q x and z = 0 Bottom temperature distribution 

t @ z=w x and y = 0 Side wall temperature distribution 
at y = 0 and z = w glass contacts with the ceramic material. The tem- 


perature of the glass and ceramic material will be same at these points. 


Following will also apply 


(aly Ke ote 
ay-y=0 ki * ay 

k oT 
oT =a Cc Cc 
Few Ea ( az) 


The boundary at y = d is the radiating boundary. In the model it is 
assumed that the enclosing surface and the glass surface behave as two 
infinite opposed parallel planes. If h is the equivalent coefficient of 


heat transfer for the convection-radiation exchange for the volume element 


dx conduction 


below 
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Energy Input Area = dx‘] 


-k' dx (34) V5 + oe,T, 


Energy Output 


hdx (T- T.) + oe Tae. 


G 
Accumulation 
. oi 4 Acme al 
~k'dx yes + oe.T. dx-hdx(T - ie) - oegl dx = oC, Sdx a 


let 6>0, soyrd 


4 h 
Se aca ) - ar aT - i) 


is the glass to air enclosure boundary condition. F is the view factor 
which accounts for the geometric arrangement of the two surfaces. With 


the assumption above following equation applies: 


] 


F => ———_ + _— 
Wee “ T/e -] 


Where Eq and Eq are the emissivities of the enclosing surface and 
glass, respectively. To go one step ahead, this definition of view 
factor is replaced with a new one, which also accounts for the re- 
radiation from the connecting walls. Glass and ceramic surface being 


gray surfaces connected by nonconducting, reradiating walls, the new 


Shape factor is given by Chapman [Ref. 3] as 





a 
1 1] 
eet} -aee- 1) 
Fiig Ag €. ae 
A, - Ay Fe 
F 1-2 
1-2 °° K¥ES- 7A, 


F is the geometric shape factor, A, and A, are the areas of the 


IZ 
radiating surfaces. 

C. NUMERICAL METHOD USED TO SOLVE THE HEAT EQUATION 

Partial differential equation describing the steady-state temperature 


T(x,y,Z) in a moving rectangular slab of molten glass in a forehearth is 


2 
k' aT Ke ee T 3 - 
Reg =o aE ae ae 3X (VT) = 0 (2 mb) 


p ay p 32 





Where cartesian coordinates are used,x is the flow direction down 
the length of the slab from the front of the hearth, y is the direction 
from the bottom channel toward the top of the slab, and z is the direction 
from the channel center toward the side. k', k", C, are the known 
constants. The initial temperature distribution at x = 0 is assumed 
to be linear in y and z. It can be found by linear interpolation from 
the four corner temperatures. Similarly it is assumed that the temper- 
ature on the boundaries (y = 0), bottom of the channel, and (z = w), 


Side of the channel, are also linear. I.e., 


T(O,y.Z) = $4(y.2) (2.2) 
T(xX,0,Z) = $5(X.2Z) (2.3) 
T(XsyoW) = $3(Xsy) (2.4) 


and are all known. The temperature of the glass is symmetric with respect 
to the center line (z = 0); therefore, 


gs 0 25) 


Radiative boundary at the surface of the glass (y = d) is that of 


radiative and convective heat transfer. Equation derived previously was, 


(2) ig = F(T,’ - 1) - Be (T- 1,) 
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Using linearizing approximation on the nonlinear part, and changing 


Sign on the second term of right hand side, 


(PE) ag = “Br (T, + 460)" (T, - Tyag) + Br (Ty = Tyeg) (2-6) 

Where o, F, k', h are Known constants. Us is the air temperature on the 
glass surface. Ty the temperature of surroundings is determined by 
linear interpolation between the temperature at the center of the channel, 
T., and the known temperature on the boundary at the side (z = w). The 
center temperature is assumed to decrease linearly from initial temper- 
ature Te at x = 0 by an amount ie: at the end of the channel (x = XL); 


therefore 


With these, the problem is formulated as being the solution of the 
partial differential equation (1), with initial condition (2.2), and 
Bewadary Comentions (2.3), (2.4), (2.5), and (2.6). 

An unconditionally-stable alternating-direction method [Ref. 4] was 
applied after passing equation (1) into finite difference form. This 
is accomplished by first getting the first derivatives with backward 
differences then substituting these in forward difference equation to 
get finite difference equation for second space derivatives, in y and 
* directions. Forward difference form is used for the derivative in x 
direction. 

The first increment step is taken as implicit in y direction, the 
second implicit in z direction, and so on, x step size being the same in 
each case. Temperatures at successive planes in x direction are related 
to each other since the first step involves values at x + Ax and these 

+] 


* 
intermediate values (Te jare used in the following step during the 


numeric calculation. System of equations obtained are, 


2] 


ant ant + 
baaTengy ~ TH) = SUTHF Lg 2TH + Tata) + SOT 91-279 aT oe) 


eps ia i-1.d Aes el 
| *n+] 
brew ee PAL Clete! wees cael pint pot nt), -nt] 
De ag Tag t Tyaq gf SOT, pet ee 


With following definitions using increments Ax, Ay, Az 


b,j = V(x Jaz) 


c. = 1/2x — 
i ells 
c. = /2ax ki 
2 nee oC 


i j= T(nax, iay, Jaz) 


Ts = T(nax, idy, Jaz) 


The indice i runs from i = 1 to i = I where i« Ay = d, indice j runs 


from j = 0 to j = J-] where jedz = w. The initial temperature gives 


sJ 
(max, Jaz) when 1= loandeile aoe. = ¢.(nax, iayiwhened) ie 
a : ce eles a2 \ a. wae 


* 
the values it az o,(y5Z); boundaries give the values Ti wes ie = 


Also, symmetry around the center line of the channel gives rise to 
and ge - lias when j = 0. Radiating boundary results 
in substitutions for That og and yom in terms of TF and ie when 
1 = [. b15 bo and 2 are obtained from corner temperatures by linear 
interpolation. 

The system of equations is solved starting from ae getting the 
intermediate solution from the first and using it in the second, making 
use of Thomas alana [Ref. 5] to solve resulting tridiagonal matrices, 


and proceeding through the channel at every y, z plane separated by 


Le 


AX in x direction, until the end of the test section is reached (x = XL). 
For integration purposes x length of the channel is divided into 500 
increments (Ax = XL/500), y coordinate is divided into 30 increments, 

and z coordinate is divided into 75 increments, making use of the sym- 


metry around the center line. 
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IT]. VeROCKIY ReRIEES 


A. PARABOLIC 

First considerations of velocity distribution within the flowing 
mass of molten glass by Duffin [Ref. 1] had involved moving side walls 
of the channel to infinity which gave velocity in X direction as a 
simple function of y coordinate. Next step was evaluation of finite 
velocity step surfaces based on a given mass flow rate of glass in an 
effort to account for higher velocities in the center channel and lower 
velocities near the side walls and bottom. None proved satisfactory, 
so final choice was the use of explicit expressions for parabolic 


velocity profiles. Physical picture is as follows: 





Bn a ll 
Zz we 


Maximum velocity occurs at y = d and on center line z = 0. Derivation 
of the general equation of parabola families gives velocity in x direction 


as a function of y and z, and physical parameters w, d and Mey as follows: 


V(¥s2) = Vay 1- (4° - (Be + 2? | (3.1) 


from which 9/4 Via, = Vege can be obtained. After integration over the 


channel cross section and dividing by channel area. 


B. VELOCITY PROFILES USING OPEN CHANNEL FLOW EQUATION 
Above parabolic profile does not account for the viscosity-temperature 


relationship. Cooper [Ref. 6] pointed out that the effects of moderate 
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viscosity gradients on velocity are negligible. Temperature in the 
forehearth changes from 2400°F to 2000°F and across this temperature 
range viscosity varies by a factor of ten. It has been shown that a 
viscosity change of 5 x 10° would affect velocity only by factor of 
three. These facts justify the assumption of constant velocity profiles 
through the channel length. 

To investigate the possible outcome when above assumption is relaxed, 
parabolic profiles are replaced with the solutions of the differential 
equation describing laminar flow in open channels. The differential 


equation to be solved is 


a 7, ou Cm eechys dpe 
ay ( ay | Seas sellee aouemmnll ae 


Where x is the direction of the flow, u is the velocity, and yu is 
the viscosity. Channel width is taken as w and depth being h. A co- 
ordinate system is set at the midstream center with boundary conditions: 

u=Q0 @ z=ew (3.3) 
u=0O0 @ yeh 


als zi 
w~=0 @ z=0 
If viscosity can be considered to be independent of depth and width, 


above equation simplifies into 


2 2 
aye oy lt @P = 0 fora 
Wile 3° aes 
which is solved with its boundary conditions by the use of infinite series. 


Solution taken from Timoshenko and Goodier [Ref. 7] is 
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Integrating over the area of half-section and dividing by the area (h-w) 


gives the mean velocity U. 


oo 
2 
= Pt ] 2h niyw 
n=1,3,5 


Data for the average mass flow rate gives an estimate for u value, 

which in turn can be used to get a value for \. An experimental value 

of u from GCIRC report [Ref. 8] was 7.77 ft/hour for green glass. This 
is used in a computer program to evaluate A taking ten terms of the 
series (3.6). Result was 0.03649. With that value, velocity at every 
node of the initial plane is calculated and held Cont through ‘the 
channel, as a first try to change previously used profiles. To give an 
idea about the difference obtained in velocity profiles, isovels are 
plotted by 3.0 ft/hour increments in Figure (3b). These changes of main 
computer program increased the running time considerably. Steady-state 
results are reported in Appendix (B), under run (A) for u = 7.77, run (B) 
for u = 8.54 and run (C) for u = 7.00. Purpose of the last two runs was 
to get an idea about the effect of mass flow rate on the temperature 
distribution, since the reported value is just a close estimate of the 
actual mass flow rate in the forehearth. Generally better correlation is 
observed to actual data in the center line temperatures, but rather poor 
results observed toward the walls of the channel. Maximum deviation of 


15 degrees F. is calculated in these three runs. 


C. VELOCITY PROFILES WITH TEMPERATURE DEPENDENCE 
To build in the viscosity-temperature dependence, available 
viscosity data, Appendix (A), is fitted to a functional form as, 


— _atbT 
=e 


This functional form for viscosity is used to get the values of u at the 
nodes for the numeric calculation and at the same time employing an 
estimated value for (dp/dx). Isovels depending on initial temperature 
distribution are shown in Figure (3c). 

To solve the partial differential equations i.e., heat equation and 
flow equation, the following route is taken: Start with initial tempera- 
ture distribution which is obtained by the use of linear interpolation 
of corner temperatures and calculate velocities at the nodes. Based on 
these velocities, temperatures in the next plane in flow direction are 
calculated by the use of the computer program already developed. All 
the physical parameters are taken as they were in the previous parabolic 
constant velocity profile computations. Computed temperatures at the 
nodes are used to recompute the velocity distribution, which is compared 
with the first velocity array until they agree within the small value e. 
Iterative technique is used in the following way: the first velocity 
array is stored in BA (i,j) array for comparison purposes; after the 
computation of temperatures in one step down the flow direction, they 
are used to get the new velocity array BIJ (i,j). If those two do not 
agree when compared at every node of the plane, a linear combination of 


the two arrays is taken as 7 
BIJ(i,j) = W.-BACi-J) + (1 - W.)BId(i.J) 


where We is a weight factor. The BIlJ(i,j) array is stored in BA(i,j) 


(| 





A. Isove's using Parabolic Equation. 





C. Isovels with Temperature Dependence 


Figure 3. Isovels for Channel Flow 
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as the new velocity array, and above cycle is repeated until convergence 
is obtained. Flow diagram is shown in Appendix (D). Numerous runs for 
the first two planes proved that the above method is appropriate. Value 
of We is varied to have an idea of its influence on the convergence and 
also on time to run. Final decision was W. = 0.25 since it cuts the 
running time by half when compared to We = 0.75. 

Due to the inevitable increase in the computer time with the above 
scheme, instead of updating the velocity array at every A x increment 
this is done at every 50 A x increments. This is updating the velocity 
array ten times through the channel length. By the trial runs to assign 
a value to c, it is found that it takes more than 18 iterations for the 
convergence with « = 0.001. Five iterations are needed when e = 0.0] 
and 216 seconds to calculate the temperatures at x = Ax and xX = 2AX. 
Taking these facts into consideration and also considering the reliability 
of the data, choice of « = 0.1] is considered appropriate. For run (D) 
velocity array is corrected with 50 A x increments, e« value was 0.1, 
running time 1/4 hour. Results at the nodes where data was taken was | 
comparable to the best run obtained with the constant velocity calculations 
but not any better. In subsequent runs (D) and (E) the estimate of (dp/dx) 
1s increased by 10% and 15% respectively, giving closer fit to the actual 
data. The absolute average deviation was around 4.6 degrees F. in all 
above runs. The largest deviation observed was at y = 5.88 inches and 
z = 8.63 inches with values around -15 degrees F. This deviation was 
also present in the results of the previous model. Another variation 
tried to compute velocity array using 25 A x increments. Results are 
reported under run (F). No improvement is observed but running time is 


increased to twenty minutes. 
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In run (G), boundary temperature at x = XL, y = d and z = w was 
decreased from 2078°F to 2035°F, and velocity updated with 50 A x 
increments. Thus run gave the smallest absolute average deviation 
4.2°F between data point temperatures and predicted temperatures. All 
above runs compute velocity array at the nodes using 10 terms of the 
defining series equation (3.5). When 5 terms of the series is employed 
running time decreased more than three minutes. The last run in these 
lines was run (H) with 5 terms for velocity calculation, « value of 0.1, 
and radiating boundary temperature at the channel corner is increased to 
2140°F from originally used 2125°F. It took 11.8 minutes to run, maximum 
deviation was 12.1°F being 1.8°F less than the best result obtained with 
the previous model. Absolute average deviation from data was 4.22°F. In 
run (I), previously used view factor is replaced with the new one, 
accounting for the reradiation from the channel walls. Maximum deviation 
from the data was 10.6°F, and absolute average deviation was 4.7°F. The 
predicted temperatures at data points for all above runs are listed in 


Appendix (B). 
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VI. PHYSICAL PROPERTY-PARAMETER ESTIMATION 


Physical properties used in the equations must be estimated or 
calculated. Many of these properties and parameters are set at values 
used in the previous model. Since this work is based on the green glass 
only, important parameters that were used in original model for that 


kind of glass are listed below. 


Radiation Conductivity Ker Ono BUyRRearT. & 
Density p 147.77 bBo, Fl. 
Specific Heat a 02375  BiURIEB. F 
Radiant Surface Emissivity Ee 0.9 

Glass Surface Emissivity Ee 0.9 

Convective Heat Coefficient h 20.0 CURR JET. iF 


In addition to the above, parameters a and b in the functional form 
ea tbT 


are used to get viscosity at the nodes as a function of tem- 
perature. Calculated values for different kinds of glass are given below. 


Plotted viscosity-temperature curve for green glass is presented in Figure 


(4). 


a b 
Green 9.07 2.82 107° 
Flint 8.60 2.60 107° 


To estimate (dp/dx) value; first A is calculated from the average velocity 
data based on the mass flow rate, then this value is used together with 


the average viscosity, since » is defined as 


hs —\ (dp/dx) 
Calculated value for (dp/dx) is 6.849 107". 


Si 


S» 
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GREEN GLASS 


Figure 4. Viscosity-Temp Curve. 
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In view factor calculations, which account for the reradiation 
from the channel walls, areas of the radiating surfaces are assumed to 
be equal. Emissivity values are not changed. Geometric view factor 


is taken from Chapman [Ref. 3] as 0.56. 


Si 


Ve DISCUSSION 


Comparison of the results obtained by the previous and modified 
models made it clear that the assumption of constant velocity distri- 
bution through the channel does not have significant effect on the 
final results. Predicted temperatures at the nodes were still off by 
small values. Model generally proved flexible and stable for all 
variations made of velocities and of boundary conditions. It still 
does not account for the possible viscosity changes due to pressure 
changes, since pressure changes with position. However, this is not 
considered significant since putting in viscosity-temperature depen- 
dence proved that initial velocities at the nodes decrease by 20-60% 
depending on the position. Further work could involve putting in 
functional forms for density, heat-capacity and thermal-conductivity 
to build in the temperature dependence of these parameters. If this 
is done, the model will be further complicated and limit its possible 
use in a closed loop computer control system. To serve for control 
purposes, time consumption should be decreased considerably. This will 
be possible only by the use of some simplifications. If the accuracy 
in positioning of thermocouples and reliability of the temperature data 
is considered, the present form of the model gives rather good agree- 
ment for the steady-state temperature data. It can possibly be used for 
design purposes. Requirement for a powerful computer is also a draw- 
back for control purposes. An alternate route might be use of a hybrid 
computer instead of a digital one to make use of the fast integration 


abilities of analog computers to decrease the running time. 
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Temperature °F 


2330. 
Ee 
Ze00'. 
Aas || 
2185. 
2155. 
cm. 
2089. 


APPENDIX (A) 


GREEN GLASS VISCOSITY DATA 
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9 
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Viscosity (Centipoises) 
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APPENDIX (C) 


FORTRAN VARIABLE DEFINITION 

N Maximum value of X direction index n 

I Maximum value of Y direction index i 

J Maximum value of Z direction index j 

XL Length of channel in X direction 

YL Length of channel in Y direction 

ZL Length of channel in Z direction 

XEND Value of X at which to terminate integration 

RHO Density of the glass 

cP Specific heat of the glass 

TKP Glass equivalent thermal conductivity in 
Y direction 

TKPP Glass equivalent thermal conductivity in 
Z direction 

SIGMA Stefan-Boltzmann Constant 

H Convective heat-transfer coefficient 

ES Emissivity of ceramic radiating surface 

EG Emissivity of the glass 

TSI Temperature of radiating boundary @ X = 0 

TO00 

sie Initial plane corner temperatures 

TODW 

TLDW Corner temperatures @ x = XL 

TLOW 

TA Temperature of the ambient gas 

TC Temperature of radiating surface along 
channel centerline 

VMAX Max. glass velocity for parabolic 
distribution 

TODWR Radiating surface Tem. @ x = 0, y= YL, 
Z See 

TLDWR ee aa surface Tem. @ x = XL, y = YL, 
Z ame 

IRNT Counter for steady-state program printout 
contro] 

NTERM Number of terms to used in defining series 
for velocity calculation 

ITR Counter to recalculate the velocity array 

ITRA Counter for number of iterations for 
velocity convergence 

ITRB Counter used for non-converging velocities 
at the nodes 

NUP Counter to update the velocity array 

WA Weight factor for linear combination 

TG 2-Dimensional temperature array 
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Bid. 47575) 2-Dimensional array for glass velocity 
(V) in X direction at (i,j) points 
in’y,z plane 


BA (k,j) 2-Dimensional array for intermediate 
velocities 
TASS. 3) 2-Dimensional array to store tem. @ the 


former plane 


se 


APPENDIX (D) 
FLOW GRAPH TO CALCULATE VELOCITIES 


ITR = 0 (Set counters to zero) 
ITRA = 0 









Z2 Calculate BIJ(K,L) 


LF 
ITR = 0 


BA(K,L) = BIJ(K,L) 





Get linear combi-| 
nation of Bld &BA 


TAS(K,L)=T(K,L) 


TITRA = ITRA + I ) 


RITE TEMPERATURE | 


= en reyes ee 


| ITR = ITR +1 | | 

nn | 

| Calculate Temp. | 

falculate Temp. : 
X= xX + AX | 

i wi 4 


| GO T0 22 


aed 





T(K,L) = TAS(K,L) 
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C COMBINED STEADY-UNSTEADY STATE GLASS PROGRAMS. 
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C READ ARRAY AND SWITCHING( PATH) DATA 
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C INITIALIZATION AND CALCULATION OF VARIOUS SYSTEM CONSTANTS AND PARAMS 
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SUME INITIALIZATION 
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CALCULATE VELOCITY DISTRIBUTION DEPENDING ON TEMPERATURE DISTIBUTUON 


0**(9.072-0-00282*T(KyL))) 
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C CALCULATE DEN2 
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ITS END VALUE 


25 IF (X-XEND) 70,105,105 
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TYPICAL DATA SET FOR STEADY-STATE CALCULATION 
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